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FOREWORD 


The work reported here concerning formulation and solution 
of a two phase mixed Stephan problem is relevant to alloy 
solidification and crystal growth processes such as those being 
investigated in low gravity experiments aboard orbiting 
laboratories. The work was supported by the Atmospheric Science 
Division of the Systems Dynamics Laboratory at the Marshall Space 
Flight Center. 
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CHAPTER I 


INTRODUCTION 

Heat transfer problems dealing with the melting or freezing of 
a substance have received a lot of attention for more than a century. 

The problem is characterized by a moving surface of separation between 
the two phases, with release or absorption of latent heat at this inter- 
face. Of interest are the position of this interface with respect to 
time and the temperature distributions in both phases. This class of 
problems are generally referred to as the Stefan problem. 

The simplest problem in this area is the one-dimensional solid- 
ification in a semi-infinite region. The first published work is by 
G. Lame and B. P. Clapeyron (1831) [17],^ which dealt with determining 
the thickness of the solid crust generated by the cooling of a liquid 
initially at constant temperature. They found the position of the 
interface to move proportionally to the square root of time, but they 
did not find the constant of proportionality. The first known solu- 
tions for this problem were found by Neumann (1860) and Stefan (1889) 

[6]. The solution takes the form of error integrals and the position 

of the interface was found to be S(t) = 2 A 4c t , where k is the thermal 

s s 

diffusivity, and A is determined by the conservation of energy equation 
at the interface. 

Further complications arise when the material freezing is an 
alloy. The classical Stefan problem normally has three unknowns, 

^Numbers in brackets refer to similarly numbered references in 
the List of References. 
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namely, the interface position and the temperatures of both phases. 

The melting temperature of an alloy is strongly dependant on the com- 
position at the interface. Therefore, freezing (or melting) can occur 
over a range of temperature. That, and the concentration distributions 
of the two phases, adds three more unknowns to the problem. Work done 
by Tao [18] considers this problem with an arbitrarily prescribed heat 
flux at the lower boundary. His closed form solution is comprised of 
an infinite sum of error functions. 

When considering a slab of finite depth, the similarity trans- 
formation, in general, cannot be applied and no known closed form 
solution exists. Boley [3] has formulated a solution for the system 
as it first undergoes solidification. However, once the effects of 
freezing have reached the upper surface the solution is no longer 
valid. To describe the system for all time, we must resort to numer- 
ical approximation techniques. Meyer [11,12] has developed an 
algorithm using the Implicit Imbedding technique that solves the prob- 
lem of the slab, but does not take into consideration the case of alloy 
solidification. 

The study presented here considers alloy solidification in a 
finite slab, with heat dissipation from both the upper and lower sur- 
faces by convection. Though the governing equations are linear, these 
interfacial conditions are not, making this a non-linear problem. 
Another complication will be added, that of allowing for density dif- 
ferences between solid and liquid phases. This causes a moving upper 
boundary as well as the moving interface. We see this solution as an 
extension of Meyer's method. Two methods of solutions will be 
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presented. One, a purely numerical technique, makes use of a Runge- 
Kutta integrator to solve the equations generated by Meyer's method. 
This solution failed at the low diffusion rates of solids because of 
the stiffness of the equations. The second solution is an analytical 
method. Closed form solutions were found of the basic equations, 
though quadrature techniques were still required. This enabled us to 
reach lower diffusion rates, though not as low as desired, because of 
accuracy limitations of the quadrature technique. 

Stefan type problems have significant importance to many indus- 
trial processes. In the production of steel billets [14], carbon 
within the metal oxidizes at the surface setting up strong concentra- 
tion gradients that cause more carbon to diffuse to the surface. This 
outer layer is no longer of the same type steel and would have to be 
ground off. Determining and minimizing the depth of this layer is 
necessary for the reduction of costs. Corrosion problems are similar 
in behavior. In the manufacture of glass [14], heat transfer rates 
must be optimized to limit crystal growth and reduce the formation of 
gas bubbles. Recently, much interest has been generated in the 
possibility of processing materials in space. In the reduced gravity 
environment, it would be possible to produce large crystals of uniform 
properties and to manufacture materials with unique properties. 
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CHAPTER II 


PROBLEM FORMULATION 


In this chapter we will discuss the formulation of the basic 
equations for the two phase mixed Stefan problem. Exact solutions by 
Neumann and Stefan and the numerical technique by Meyer will be 
reviewed. For the mixed Stefan problem we will consider a slab con- 
sisting of a binary liquid alloy, initially at constant temperature, 

T q . The slab is of finite depth, d^, and infinite width (Figure 1). 
Heat is to be removed from both the upper and lower surfaces, with 
solidification starting at the lower surface. The problem is charac- 
terized by a moving interface, S(t), between the solid and liquid 
phases, with latent heat release at the interface. The equations will 
allow for different densities for the solid and liquid causing a 
volume change and therefore a moving upper surface during freezing. 

The upper and lower surfaces will be closed to any mass flow and mass 
diffusion in the solid and liquid will be investigated. Since the sub- 
stance is an alloy, the freezing temperature will be dependant on the 
composition at the interface. A sketch of this model is shown in 
Figure 1 . 

The governing equations for the temperature and concentration 
distributions when there is a change in volume during solidification 
(p /p ^ 1) take the form [6]s 

S x, 


9T„ P s dS 3T 
3t dt 3z 


= < 


3 2 T, 


3z 
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Figure 1. Sketch of theoretical model. 
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in the solid. C denotes the solute concentration and T the temperature. 

S (t) is the interface position and p, k and d denote the density, thermal 
diffusivity and solute diffusivity. The subscripts s and Jt, denote solid 
and liquid, respectively. 

The interface conditions are: 


k 

s 



_ . fh ds 

i 3z P s^ dt 


(5) 


P d 
s s 


^8 

3z 


- P A 


3z 


- p s C * (1 


k) ^ 
k) dt 


( 6 ) 


Equations (5) and (6) are the conservation of energy and mass equations 

at the interface, k, k and l are the thermal conductivity, the partition 

coefficient and the latent heat. Also at the interface, both the solid 

and liquid must be at the equilibrium melting temperature, T , and the 

m 

proper value of solute concentration which may be described by 
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( 7 ) 


T = T, - mC ft 
mi H 


C = kC 0 
s l 


( 8 ) 


Equation (7) implies that the melting temperature of the alloy 
is a linear function of concentration. Looking at the phase diagram. 
Figure 2 [1], we can see that this is a simplification of the actual 
relationship. The decrease in the melting temperature with increase 
in concentration can be rather closely approximated by a linear func- 
tion. However, at a certain concentration the melting temperature 
will begin to increase. Since the method of solution does not depend 
on this equation, a more complicated expression could be substituted, 
but will not be for this work. 

The boundary conditions at the upper and lower surfaces take 
the form: 

at z = 0 


k s f = h l < T s - V 


(9) 


3C 

3z 


( 10 ) 


at z = d 


0 



- 1 ) S(t) 


-k 

«, 3z 


- (T. 


- T 2> 


( 11 ) 
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Figure 2. Pb-Sn Phase diagram. 
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8C 
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where subscripts 1 and 2 denote lower and upper surfaces, respectively. 
Equations (9) and (11) represent heat dissipation from the surfaces by 
convection into a medium at temperatures and T^, respectively. 


The Solutions of Neumann and Stefan 

The first solutions for the Stefan problem, in a semi-infinite 

space, were found by Neumann in 1860 and Stefan in 1889 (6]. The 

governing equations are (1), (2), and interface condition (5). They do 

not allow for a density change between the solid and liquid, i.e., 

p = p . The boundary conditions at zero and are 
s X- 

T = T at z + ® 

T = 0 at z = 0 
s 

The solutions 


T = T - Berfc ( — - — ) 

* 2/Tt 

s 


T = Aerf ( — - — ) 

S 2/TT 

s 


satisfy the basic equations and boundary conditions at zero and <», 
where A and B are constants of integration. Using the interface con- 
dition 
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we have 


T =T=T at z = s 
s l m 


Aerf ( — - — ) = T - Berfc ( — - — ) 


= T 


2^t 


2/k t 
s 


m 


For this to be valid for all t, we must have 


= A = constant 


2/k t 
s 


or 


s - 2A/TT 
s 


rewriting the solution at the interface, we have 


Aerf (A) - T - Berfc (A 


j%. 

J K„ 


m 


which determines the constants A and B as follows: 


m 


erf (A) 


T - T 


m 


erfc (A 


v k/ 


Now, using the interface condition (5) , we can evaluate the constant X 
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exp(- 


) 


exp(-X 2 ) 
erf (X) 



k 

s 



T - T 

m 


T 

m 


‘s X 


erfc(X 



Xfc/ir 

C T 
P s m 


Stefan's solution assumes the liquid is initially at the melting 
temperature (T = T m ) , giving the equation for X as 

C T 
P s m 

Xexp(X )erf(X) = 

ji/ir 


Adding the change in density term (p /p - 1) » so that the basic 

S X# 

equation matches equation (1) , changes the solution in the liquid 
region to the following: 


T 


s 


= T - Berfc 


— + X(— - 1) 
2»^t p £ 



with 


B 


T 


- T 

m 


erfc(X- 



and the equation to determine X: 
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(T - T )k exp(- 

m Jo 


exp(-A^) 
erf (A) 


T k erfc(A 
m s 


.2 2 
A p K 
s s 


p * 


J? 


fs V 

k a p * 



X&/jT 

c T 
p s m 


Numerical Technique, the Method of Meyer 

When considering a slab of finite depth, a numerical technique 
must be employed. In a series of papers Meyer [11,12] has developed 
one such technique using the Invariant Imbedding Method. This tech- 
nique also allows for the use of non-linear boundary conditions. The 
governing equations are (1) and (3), with interface conditions (5), 
and boundary conditions (9) and (11). Since numerical methods are 
used, it is convenient to first nondimensionalize the equations. The 
following substitutions were made: 


R 


C 



D 


d 

s 



Lewis Number 


M 



P 

s 


d 

— = Lewis Number 

K 


K 



Vo 


= Biot Modulus 
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z = 


L = 


Cp < T 1 - T o> 

JO 


A - 17- 


Y = 


T 2 ~ T 0 
t i - T o 


0 = 


T - T 
_ ^0 

T - T 
1 0 


T 

d o 


After nondimensionalization, the governing equations, interface 
and boundary conditions take the form 


90£ - (R - 1) ^ - 9 ° Z 


3x 


dt 3z 


3z 


(13) 


3T 


3 2 0 


= K 


3z 


(14) 




(15) 


30 (0) B. 

-Sr- - t ». - » 


(16) 


30 £ (Z f ) 

-IT- - - B 2< e * + » 


(17) 


Meyer’s method consists of first discretizing the time operator 
in the following manner: 
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3t At 

yielding a set of ordinary differential equations, solvable at discrete 
time steps, N. Next, equations (13) and (14) are rewritten as a first 
order system. 


U 


£ 


= 0 


£ 


^ N 

9 £ - 6 * 
At 


N-l 


N N-l 
(R - 1) (S a - S L J 

At 


U n = U' 
£ £ 


U = 0' 
s s 


(18) 


0 


N 


- 0 


N-l 


At 


= KU' 


(19) 


where prime denotes differentiation with respect to z. With the con- 
servation of energy equation at the interface taking the form 


A 




RL 


s" - s 11 - 1 

At 


( 20 ) 


Now, the Invariant Imbedding technique is used to solve for 

and U . The solution of D. is assumed to take the form 
s £ 


N N N N 

U” = Y ©„ + Z„ 

£ £ £ £ 


then 
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N' N' N N' N’ 

U. = Y + Y 0' + Z “ 

£ £ £ £ £ 


N' N' N N N' 

U„ = Y. + Y„ U + Z„ 

£ £ £ £ £ 


These equations are substituted into (18) to yield 


YfV + (Y* 1 ) 2 0„ N + y' , Z , N + Z, N ’ 
£ £ £ £ £ £ £ 


N 


Ax 


0 


N-l 

£ N N N 

JT - V* s ? - V* 


where 


(R - 1) (S N - S N_1 ) 
F R Ax 


rearranging, 


y *\ n + + f r y z < ■ 


- z, N ' - y, n z, n - f„z, n -^- 

£ £ £ R £ Ax 


from this we get the equations for Y and Z . 

JO JO 


N ' = _ rv N i 2 - F y n + -L 

£ (Y £ } F R Y £ + Ax 


( 21 ) 


N' 


N-l 

_ N N N _£ 

" “ Y £ Z £ “ F R Z £ ' Ax 


( 22 ) 


The boundary conditions are derived similarly. 


N N N N 

U = Y„ 0„ + Z” 

£ £ £ £ 


is substituted into (17) to yield 
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(Z f ) 0 £ n (z f ) + z £ n (Z f ) = - b 2 (0 £ n (Z f ) + Y) 

and therefore 

Y * B <Z f> - - B 2 

Z l < Z £> ' - V 

where = 1 - (R - 1)S^. Similarly, 


yields the following equations and boundary conditions: 


dY 


N 


dz 


JL_ ( y n ) 2 

KAT K 8 ’ 


(23) 


at z = 0 Y 


N 


dZ 


N 


dz 


- y n z n 

s s 


N-l 


KAt 


(24) 


at z = 0 Z 


N 


After solving equations (21)-(24), the temperature distributions can 
be obtained by solving the following equations: 
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(25) 



dz 




+ Z 


N 

i 



dz 



N 


+ Z 


N 


(26) 


Meyer's Method for the Mixed Problem 

Meyer's method can easily be applied to alloy solidification. 
Equations (2) and (4), interface conditions (6), (7) and (8), and 
boundary conditions (10) and (12) are nondimensionalized, time operator 
discretized, and then written as a first order system yielding: 


V„ = C'„ 
1 l 


C N - C N_1 

JL 

At 


(R - 1) (S N - S N-1 ) 
At 


V = P V ' 
l it 


(27) 


V = C' 
s s 


c N - c 1 *- 1 

S S 

At 


KP V f 
s s 


(28) 


dC dC 


RD 


dz 


.foaL - . s*- 1 

dz P <! 1 At 


0 = 0. - MC„ 

mi l 


C = kC n 
s S, 


(29) 

(30) 

(31) 
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( 32 ) 


V and V are assumed to have solutions of the form 
a s 


N N N N 
V * = X * C * + R * 


N N N N 

V » X C + R 
s s s s 


N N N N 

Then the equations to solve for X^ , X g , and R g are 


F R x n 
* 


(xf ) 2 


at z = Z f X" 


r^ + x» 1R » 


at z = Z f R " = 0 


dX , „ o 

s_ _ 1 _ ,„N. 2 

dz " KP At V s ' 
s 


at z = 0 X = 0 
s 
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(37) 


dR 


dz 


- - X 


N r n 


s 


s 


,15-1 


KP AT 

s 


N 

at z = 0 R = 0 
s 


After solving equations (34) — (37) » the concentration distributions can 
be obtained by solving the following equations: 



dz 





(38) 


dC 


N 


s 

dz 



+ R 


N 


(39) 


Calculation of Initial Condition 

In order to solve equations (22), (24), (35) and (37), temper- 
ature and concentration distributions from the previous time step are 
needed. This necessitates initial distributions for both the temper- 
ature and concentration at the onset of solidification. For the present 
case, we will assume that no mass diffusion occurs prior to freezing. 
Thus, the concentration distribution is taken to be 

C £ (z, t) = 1 

throughout the liquid, up until solidification starts. 

Since the initial temperature, 0^, may not initially be equal to 

the melting temperature, 0 , such an assumption may not be made. We 

m 

must therefore derive an equation for the temperature as a function of 
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position and time, which is appropriate to the specific cooling at the 
upper and lower surfaces. The governing equation [6] for this is 




(40) 


with boundary conditions 

IT" V s * - D ' 0 < 41 > 

9e t 

“aF “ ■ V e * + ,) at z = 1 (42) 

and with initial conditions 0 = 0 at t = 0. The solution is derived 

■JO 

using separation of variables. The solution of the homogenous part 
takes the form 


Q z (z, t) = Z(z)T(t) 


then 


ZT ’ = Z"T 


or 


T 


zl' 

z 


- X 


2 
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f - -x 2 t 

dx 


^ - -X 2 Z 
dz 


T - A 1 exp(-X t) 


Z = A 2 sin(Xz) + A^cos(Xz) 


and 


0^ = (^sin^z) + A^cosCXz)] exp(-X t) 


(43) 


30 


Z 2 

= X[A2Cos(Xz) - A^sinCXz)] exp(-X t) 


(44) 


at z = 0 


90 £ 

~dT " V* 


2 2 
XA 2 exp(-X t) = B^A^exp(-X t) 


B 1 

A 2 X A 3 


(45) 


at z = 1 


90 fc 

IT " - V* 


X[A 2 cos(X) -A^sin^)] exp(-X^f) = - B2tA2sin(X) i-A^cos^)] exp(-X^x) 


using (45) 
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sin(X) + A^cosCX) ] 


B 1 A 3 


BjA,jCos(X) - XA^sin(X) - - B 2 [ — - 


X - b.b 2 

(B 1 + B 2 )cos(X) = [ sin(X) ] 


X N (B 1 + B 2> 
tanCX ) = - 

X N" B 1 B 2 


( 46 ) 


So the solution of the homogenous part is 


0^(z, x) = Z A^ [ y- sin(X N z) + cos(X N z) ] exp(-X^x) 


N=1 


N 


(47) 


Where the eigenvalues, are determined from (46). 

For the non-homogenous part, assume the solution takes the form 


0 = a + 3z 

P 


Then, using (41) and (42) we get 


3 - B L (a - 1) 


and 


3 = - B 2 (a + 3 + V) 
combining and solving for a. 
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Bjd + B 2 ) - 
a = Bjd + B 2 ) + B 2 

Now for the total solution we have 


0 £ (z, t) = a + £z + Z ^ [ Y~ sin(X N z) + cos(X N z) ] exp(-X^r) 

N=1 N 


Now using the initial condition 

0=0 at t = 0 

X / 

we can solve for the A.. 


"1 

0 = a + Bz + Z Ajjf — sin(X N z) + cos(X N z)] 
N-l N 


- a - gz 


Z 

N-l 


¥»<*) 




/q (-a - BzjEjjjCzJdz 
'o B 2 b MdZ 




B, o aB. (aB.-B) 

-2X N { [a + Pd +-f )] + [^f- (i - Bl ) -T 7 ]eos(X N ) + -4— } 
A . T N N N 


[ Bj + X 2 n +h ( ■ ■ ~ 2 - 1 ) sin (2X n ) + 2B lS in 2 (X N )] 

X N 
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Now, to calculate the distribution for the first state, N = 0, we must 


find t = T so i suc h that the temperature at z = 0, is equal to the melt- 
ing temperature, 0^. The initial temperature distribution is then 
equal to 0(z,t p . 
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CHAPTER III 


SOLUTION AND RESULTS 
Numerical Method 

In order to solve for the temperature and concentration distri- 
butions in both the liquid and solid phases, equations (21) -(24) and 
(34)-(37) must first be solved. These equations comprise of a set of 
coupled, nonlinear ordinary differential equations with their proper 
boundary conditions. Such a set of boundary value problems are most 
suitably solved with a Runge-Kutta-Fehlberg [10] integrator. For the 

present problem, a RK7(8) was used. Before integrating these equa- 

N N 

tions, a new interface position, S , and concentration, , must be 
assumed. Once these equations are solved, the temperature and con- 
centration gradients at the interface may be obtained from equations 
(25), (26), (38) and (39). The conservation of energy (20) and mass 

(29) equations at the interface are then used to check these results. 

N N 

If the initial guess on S and was in error, these equations will 

N 

not be satisfied. A better estimate for the correct values of S and 
N 

C can be obtained from these equations through a Newt on-Rhap son 

X/ 

method. This process of evaluating the temperature and concentration 
gradients and reestimating the interface position and concentration 
is repeated until convergence to the correct values occurs. Now, 
equations (25), (26), (38) and (39) can be integrated for the tem- 
perature and concentration profiles from the interface to the lower 
surface for the solid and from the interface to the upper surface for 
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the liquid. The time is then advanced by one time step. Ax, and the 

process repeated for the new interface position. 

This method (Antar [2]) gave satisfactory results down to a 

-4 

Lewis number, P , of order 10 . At lower Lewis numbers overflow 

s 

problems occurred due to the stiffness of equation (36). This can 
easily be seen by using the following values as an example: 

P = icf 5 

s 

K = 2.7 
At = 1 

To demonstrate the problem, we will use a fourth order Runge-Kutta 
technique as outlined in ref [4], and the above values with the 
boundary condition of zero, 
if 


y' = f(z,y) with a <_ z _< b 
y(a) = boundary value 

b a 

h = — — — where N = number of intervals 


then an approximate solution to y is w, found in the following manner: 

Wg = boundary value 

k 1 = hf (x i , w ± ) 
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tl 1 

k 2 = hf (x i + j , w ± + -jkp 
k 3 = hf (x ± + y , w i + |-k 2 ) 

k 4 = hf(x i+l » w i + k 3> 

w i+l = ~6 ^ k l + 2k 2 + 2k 3 + k 4^ 


for each 1=0, . . . ,N-1. 

Applying this to equation (36) with 50 intervals we arrive at the 
following values: 

For i = 0 

k x = +7.407 x 10 2 
k 2 = - 2.003 x 10 3 
k 3 = - 1.931 x 10 4 
k 4 = - 7.460 x 10 6 

and 

w x = - 1.250 x 10 6 


For i = 1 


k = -3.127 x 10 10 
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1 ft 

k 2 = - 4.889 x 10 
k 3 = - 1.195 x 10 35 
k 4 = -2.865 x 10 68 

The calculation for causes overflow on the VAX 11/780, which has an 
upper bound of order 10^. 


Analytical Method 

Inspection of equations (21) and (34) reveals them to be of 
the following general type: 


dz 


A 


- By - y 


2 


There exists a closed form solution [16] for this type of equation. 
Finding closed form solutions avoids the use of numerical integrators, 
therefore rendering the problem of stiffness more manageable. The 
solution is found by using the transformation 

u' 

y = V 


The equations then take the form 

u" + Bu f - Au = 0 

whose solution is found in the following manner. 
The characteristic equation is 
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+ Br - A = 0 


r j 2 = | ( _B ± 2 + 4 A) 

u = C^exp(r^z) + C 2 exp(r 2 z) 


then 


where 


r^expCr^z) + Cr 2 exp(r 2 z) 
exp(r^z) + Cexp(r 2 z) 


(r x - y) 

c ° ' (r 2 - y) e * p((r l - r 2 )z) 


Similarly, equations (23) and (36) have the form 


iz = 

dz 


A - 


2 

y 


and can be solved using the same transformation 

u' 

y = V 


yielding 


u" = Au 

u = C^exp(/Az) + 
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with 


y = /X [ ex P (i/Az) - Cexp(-/Az) ^ 
exp(*/Az) + Cexp(-/Az) 


C = [ — ^ ] exp(2v^Az) 

/A + y 

Equations (22), (24-26), (35), (37-39) are of a more complicated type 

|Z = - F (z)y - G(z) 

but can be easily solved using the Variation of Parameters technique 

[9]. 


Z 2 Z ' 

y = exp(- / F(z)dz) t C - / G(z')exp ( / F(z)dz)dz' ] 

Z i Z i Z i 

Equations (22), (24), (35) and (37) contain functions whose 
values are known only at certain points, requiring quadrature evalu- 
ations to determine their solutions. These techniques, though, are 
known to be highly stable, allowing solutions to be obtained at lower 
Lewis numbers. The same iterative technique as described earlier is 
used to determine the interface position and the temperature and con- 
centration profiles at each time step. 

Now the solutions of equations (21) — (24) and (34)-(37) are: 


Y 


N 

4 


r^exp(r jz) + C r 2 exp(r 2 z) 
exp(rjz) + C 1 exp(r 2 z) 


(48) 
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r = — [ - F + /p^ + 1 

1,2 2 L R " >/ R At J 

(r l ' B 2 ) 

C 1 = " (r 2 - B 2 ) exp((r l ’ T 2> z e> 

N-l 

z ^ z 0 z 

z” = exp(-J (F R +Y £ N )dz) t C 2 -J f exp(J f (F R +Y £ N )dz)dz * ] (49) 

z z z' 


C„ = - B, 


exp(6z) - C,.exp(-Sz) 

Y w = x r - 

s L exp(<Sz) + C,.exp(-6z) 


] 


(50) 


6K - B 
] 5 = 6K + B, 


1 


6 = 


v'KAx 


Z N = exp (-/ Y N dz) [ C, - f 
0 s 6 0 


z 0 


N-l 


s 


KAx 


exp( J 


Y N dz)dz' ] 
s 


c 


6 



(51) 



r 1 exp(r 1 z) + C 3 r 2 exp(r 2 z) 
exp(r 1 z) + C 3 exp(r 2 z) 


(52) 
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r . -Lr-Is. hhyi + _*_ 

r l,2 2 Pj-J'p/ + P t ix 


C 3 ' -^ ex P <(r l - r 2 >z f ) 


/ Zf <^ + X t N )dz) ; Z f ^-exp(- j Zf (^ + x;>d Z )d 2 ' 

z & z & z & 


X N tanh(-$2-) 

s /p“ /F" 

S S 


C 

= ( 1 )/ K P At cosh ( — )dz 

cosh(— ) 0 KP s AT /T~ 

/V 3 

s 


Now the solution to equations (25) , (26) , (38) and (39) can be found 
using the Variation of Parameters technique [9]. 


0^(z)=exp(J Y^dz)[C 9 +/ Z^exp(-J Y^dz)dz' ] 


0^ (z) = exp( J Y^dz)[C 10 +J Z^exp(-J Y^dz)dz' 
z. z ± z ± 


C^(z)=exp (/ X^ ) [ C n +/ exp (-/ X^dz)dz' 

Z i Z i Z i 


C^(z)=exp(J X^ ) [ C 12 +/ exp (-/ X^dz)dz r ] 

Z i Z i Z i 
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The constants, and C^, represent the temperature and con- 
centration of the liquid at the interface. The constants, and 

C 12* re P resent the temperature and concentration of the solid at the 
interface. and must be equal. Since only the conditions at 

the interface are known, the integrations in equation (56) and (58) 
are performed from the interface up, and, in equations (57) and (59), 
from the interface down. 

Asymptotic approximations were used whenever possible to 
reduce the chance of overflow errors. For example, equation (50) is 
evaluated as 


Y 


N 

s 


= 6 


whenever 


6z > 60 


Also, exponents were combined as much as possible. In equation (59), 
for example. 



was evaluated first, then treated as a constant so that it could be 
carried within the integral and added to the exponent before the 
exponents are evaluated. 

Many of the integrals could not be solved exactly because they 
contain the concentration or temperature functions whose values are 
known only at certain points. They were therefore solved using a 15 
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point Gauss-Legendre quadrature, with maximum interval size of 0.1 
and a cubic spline interpolation [4] between the known points. 


b ,, v 15 r , (b - a) 4* b + a 

/ F(x)dx = ~ a; Z w ± F ( — ^ ) 

a i=l 


The roots and weights (r^, w^) were obtained from Table 2.2 in ref 
[5]. 

This method gave good results down to a Lewis number, P , of 
order 10 Below this value, problems occur due to the accuracy of 

the integrator, when applied to equation (59). The error occurs when 


/ 

Z i 


r n 

s 


exp (-/ X^dz)dz' 

z. 

l 


(60) 


is evaluated. To demonstrate this, using the functions 


and 


„ N 6 . . , 6z 

X = tanh( ) 

s /p" sr 

s s 


R N = 


cosh( ) 

/p~ 

s 


cosh( ) 

, z /P s 

; o 


r n = 

s 


-6 


tanh ( ) 

/?“ 
s 

/ P~ 
s 
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we evaluate equation (60) in closed form. Note that these equations 
are the same as (54) and (55) , except for the missing concentration 
term in (55) . The exact solution of (60) then is 

sech( — ^— ) - 1 

/P~ 

s 

when evaluated from zero to one. For large values of 6//P~ the solu- 

s 

tion approaches the value of -1.0. Evaluating equation (60) numer- 
ically, using a 15 point Gauss-Legendre quadrature with 10 intervals, 
and the following parameters 


K = 2.7 


P 


s 


10 


-7 


At = 1 


z = 1 


we get the value of -1.36, which is a difference of .36 from the true 

value. Table 1 lists the results of different interval sizing at 

different values of P . 

s 
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Table 1 


The Effects of Different Interval Sizing on the Accuracy 
of a Gauss-Legendre Quadrature Technique 


Value 
of P s 

Number of 
Intervals in 
Quadrature 

Difference 

from 

Actual Value 

io - 6 

10 

2.33 x 10 -2 

r^. 

1 

o 

rH 

10 

3.62 x 10" 1 

io - 7 

100 

1.22 x IO -4 

10- 7 

250 

1.30 x 10~ 7 

io - 8 

100 

3.15 x 10 -2 

io ' 8 

250 

5.24 x 10 4 

io - 8 

500 

o> 

O 

X 

f— * 

o 

1 

O' 


Reasonable accuracy is obtained for P g of order 10~ 7 using 250 
intervals. However, combined with a 15 point quadrature formula, 
this means 3750 function evaluations (actually, equation (55) has an 
additional integration, requiring 56,250 evaluations), which means a 
high probability of round-off error occurring during the computations. 
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Computation Results 


As previously discussed, the numerical technique gave good 

results down to a Lewis number, P g , of 10 . Alloys having this high 

rate of diffusion in its solid phase are rare. For more common 

alloys, such as lead-tin (Pb-Sn) , the Lewis number is of order 10 ^ 

[7]. Because of this restriction the analytical technique was 

derived. With this technique, we were able to get results down to 

—6 

a Lewis number, P , of order 10 . The following parameters, based 

s 

on the temperature differences 



were used in the calculations. 
K =2.326 
M = - 0.22 
L = - 27.7 
P t - 10- 2 
k = 0.580 
A = 2.703 



B = 1.00 
b 2 = 0.01 

R = 1.00 or 1.06 

-4 -6 

P = 10 10 

s 

¥ = - 0.14286 

At = 1 
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-4 

A more realistic liquid Lewis number, P^, is of order 10 [8]. Compu- 

-2 

tat ions were performed at a Lewis number, P^, of 10 in order to keep 

at least the ratio of liquid to solid diffusion rates realistic. Com- 

-4 

parison of the two techniques at P =10 is shown in Figures 3, 4 

s 

and 5. The interface position and the temperature and concentration 
at the interface as a function of time are shown to match excellently. 
The solid lines represent the analytical solution while the circles 
represent the numerical solution. 

The next series of plots, Figures 6, 7 and 8, show concentra- 
tion distributions at different Lewis numbers at time step intervals 
of 25, with the first time step at 25 and the last at 200. The inter- 
face position is seen as the vertical lines with the solid phase to 
the left and the liquid phase to the right. The effects of lower 

-4 

Lewis numbers are very apparent. In Figure 6, at Lewis number 10 , 

there is significant diffusion in the solid, while in Figure 8, at 
Lewis number 10 almost no diffusion is noticeable. The only 
movement is at the lower surface, which has had more time for diffu- 
sion to occur. 

Figures 9-20 show how the interface, concentration and tem- 
perature vary with time and the temperature with interface position 

-4 -5 -6 

at the three Lewis numbers, 10 , 10 and 10 . We can see in 

Figures 9, 13 and 17, that the velocity of the interface is decreas- 
ing as the alloy solidifies. The upward motion of the interface has 
all but stopped. This is due to equation (7), which determines the 
melting temperature of the alloy. As seen in the phase diagram, the 
melting temperature decreases, as concentration increases, to a 
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PS= 



Figure 4. Comparison of temperature vs. time for numerical and 
analytical techniques. Solid = analytical. Circles = numerical. 
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Figure 5, Comparison of concentration vs. time for numerical and 
analytical techniques. Solid = analytical, Circles » numerical. 
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Figure 7 . Concentration distributions at 
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Figure 8. Concentration distributions at P 
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Figure 9. Interface vs. time at P 
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Figure 10. Concentration vs. time at 
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Figure 13. Interface vs. time at P 










certain value and then begins to increase. Allowances for this have 
not been made and since the concentration of the liquid will be ever 
increasing, the melting temperature will continually decrease, which 
can be seen in Figures 10, 12, 14, 16, 18 and 20. It is also 
interesting to note that the velocity of the interface is proportional 
to the Lewis number. The lower Lewis numbers cause slower solidifica- 
tion. Due to the slower diffusion rates, the interface concentration 
is higher and therefore the melting temperature is lower, requiring 
more time for the alloy to reach the melting temperature since the 
heat rates remain the same. 

Figures 21-23 illustrate the effects of different densities 
in the solid and liquid phases. The interface moves at a slower rate 
when the density of the solid is greater than that of the liquid. 

The temperature and concentration with respect to time does not change 
significantly. 

Figures 24-29 show the overall temperature and concentration 
distributions for the different Lewis numbers. The concentration 
distributions are represented by the solid lines and the temperature 
distributions by the dashed lines. 

In Figure 30 we see the comparison of results obtained from 

the analytical method at P = 10 ^ (Figure 30-A) with intuitive 

s 

results [13] deduced from the phase diagram to the limit of zero 
dif fusivity in the solid (Figure 30-B) . 
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Figure 21. Comparison of interface vs. time for R = 1.00 and R = 1.06. 











Figure 25. Temperature and concentration 
distributions for time = 125 - 200 with P„ = 10“ 
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Figure 27. Temperature and concentration 
distributions for time = 125 - 200 with P = 10“ 
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